# Escoamento completamente desenvolvido em condutos de secao circular
#
# regime laminar (mu constante) ou turbulento (mu = mu(r))
# campo gravitacional nulo
# fluido newtoniano
# discretizacao por volumes finitos
# code: FVMpipeflow.m - octave
# para rodar: run FVMpipeflow.m
#
# Diego S. Rodrigues
# 13 de agosto de 2012
# PPG-CCMC - ICMC-USP
# Disciplina: SME-5802 - Introducao a Mecanica dos Fluidos Computacional
# Prof. Gustavo C. Buscaglia

clear all;

mu_const = 1.0;  # viscosidade dinamica do fluido
p0 = 1.0;  # pressao em z=0
pL = 0.0;  # pressao em z=L
L = 1.0;   # comprimento do tubo
r = linspace(0,1,51); # r(1) = 0.0, sempre!
w=zeros(50,1);
b=zeros(50,1);
gradp = -(pL - p0)/L; # dp/dz
nunk = length(r) - 1; # numero de intervalos (ou numero de incognitas - nunk)

for iunk=1:nunk
h(iunk) = r(iunk+1) - r(iunk);
end

for i=1:length(r)
mu(i) = mu_const;
end

for it=1:1000

gradp = -(pL - p0)/L*cos(0.1*it); # dp/dz


iunk=1;
A(iunk,iunk) =  mu(iunk+1)*r(iunk+1) / (r(iunk+1)+ 0.5*h(iunk+1));
A(iunk,iunk+1) = - A(iunk,iunk);
A(iunk,iunk) = A(iunk,iunk)+r(iunk+1)^2*10;
b(iunk) = gradp*0.5*r(iunk+1)^2+r(iunk+1)^2*10*w(iunk);

for iunk=2:(nunk-1)
A(iunk,iunk-1) = -mu(iunk)*r(iunk) / (0.5*(h(iunk) + h(iunk-1)));
A(iunk,iunk+1) = -mu(iunk+1)*r(iunk+1) / (0.5*(h(iunk+1) + h(iunk)));
A(iunk,iunk) = - A(iunk,iunk-1) - A(iunk,iunk+1)+(r(iunk+1)^2 - r(iunk)^2)*10;
b(iunk) = gradp*0.5*(r(iunk+1)^2 - r(iunk)^2)+(r(iunk+1)^2 - r(iunk)^2)*10*w(iunk);
end

iunk = nunk;
A(iunk,iunk-1) = -mu(iunk)*r(iunk) / (0.5*(h(iunk) + h(iunk-1)));
A(iunk,iunk) = - A(iunk,iunk-1) + 2*mu(iunk+1)*r(iunk+1) / (h(iunk))+(r(iunk+1)^2 - r(iunk)^2)*10;
b(iunk) = gradp*0.5*(r(iunk+1)^2 - r(iunk)^2)+(r(iunk+1)^2 - r(iunk)^2)*10*w(iunk);

#b = b'; # transformando b num vetor coluna



# Resolvendo o sistema linear
fprintf(stdout, "\n >>> Solving Linear System ... \n \n");
w = A\b;
fprintf(stdout, "\n >>> Linear System was solved ... \n \n");

# Definindo vetor auxiliar z para plotagem
iunk = 1;
z(iunk) = 0.0;

for iunk=2:nunk
z(iunk) = r(iunk) + 0.5*h(iunk);
end

#z(nunk+1) = r(nunk+1);
#w(nunk+1) = 0.0; # no slip (nas contas acima ja usado p/ iunk=nunk)

# Solucao exata
#for i=1:(nunk+1)
#wexact(i) = -gradp*(r(length(r))^2-r(i)^2) / (4.0*mu_const);
#end
#wexact = wexact';  # transformando wexact num vetor coluna

# Imprimindo saida num arquivo
#fid = fopen("FVM2.dat","w");
#fprintf(fid,"# z velw_exact velw");

#for i=1:nunk+1
#fprintf(fid,"\n %g %g %g",z(i),wexact(i),w(i));
#fflush(fid);
#end

plot(z,w,"-");
pause;

end

fclose(fid);

fprintf(stdout, "\n >>> END COMPUTATION \n \n");
